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The  only  complete,  analytic  solution  for  conduction  problems  with  phase  change  is  the  Neumann  solution.  The  Neu¬ 
mann  solution  is  valid  for  phase  change  in  a  semi-infinite,  homogeneous  medium  with  a  step  change  in  surface  tempera- 
tu restarting  from  an  initial  temperature  which  can  be  different  than  or  equal  to  the  fusion  temperature  of  the  medium. 

The  Neumann  solution,  when  applied  to  soils,  forms  the  basis  of  a  number  of  formulae  for  calculating  the  depths  of 
freezing  or  thawing..  Widely  used  graphs  were  previously  developed  that  are  valid  only  when  the  ratios  of  the  thermal 
|  conductivities  and  thermal  diffusfvlties  of  the  frozen  and  thawed  soils  are  unity.  In  this  report  general  charts,  applicable 
to  any  property  ratios,  are  developed.  The  figures  have  been  drawn  specifically  for  soil  systems,  but  they  are  applicable 
to  any  material  with  appropriate  property  ratios.  <f~- _ _ 
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PREFACE 


This  report  was  prepared  by  Dr.  Virgil  Lunardini,  Mechanical  Engineer,  Applied 
Research  Branch,  Experimental  Engineering  Division,  U.S.  Army  Cold  Regions  Re¬ 
search  and  Engineering  Laboratory. 

The  study  was  conducted  under  DA  Project  4A161 101A91D,  In-House  Laboratory 
Independent  Research. 

The  author  expresses  his  appreciation  for  excellent  reviews  by  Dr.  R.  Berg  and 
G.  Phetteplace,  CRREL. 
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NOMENCLATURE 


b  s/ITj”  =  \Af/*, 

C  volumetric  specific  heat  or  heat  capacity 

lf,  lt  air  index  for  freeze  and  thaw 
/$  surface  index 

k  thermal  conductivity 

L  volumetric  latent  heat 

n  ratio  of  surface  to  air  index 

q,  r  property  ratios  defined  in  text 

ratio  of  thermal  conductivity  of  ice  to  water 
t  time 

Tf,  Tm  r$  fusion,  initial,  and  surface  temperatures 
x  volumetric  fractions  of  soil  solids,  liquids,  gases 

X  depth  of  phase  change 

«  (Tf-TjHTt-Tf) 

y  phase  change  parameter 

0  length  of  thaw  or  freeze  season 

k  thermal  diffusivity 

X  phase  change  parameter 

(CJ2L)(Tr  Tf)  thawing 
"  (Cf/2L)(  7f-  fj)  freezing 


Subscripts 

f  frozen 

g  gas 

i  ice 

C  liquid 

s  solid 

t  thawed 

w  water 
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THE  NEUMANN  SOLUTION  APPLIED 
TO  SOIL  SYSTEMS 


V.|.  Lunardini 


INTRODUCTION 

One  of  the  most  common  problems  facing  engineers 
in  the  cold  climate  regions  of  the  world  is  the  need  to 
estimate  the  depth  of  freeze  or  thaw  of  a  soil  system. 
Despite  the  importance  of  moving  boundary  problems 
of  conduction  with  phase  change,  there  are  only  a 
handful  of  exact  solutions.  One  of  these  solutions, 
which  includes  both  latent  and  sensible  heat  of  the 
soil  mass,  is  that  of  Neumann  (ca.  I860)  which  was 
generalised  by  Carslaw  and  |acger  (I9S9). 

At  the  start  of  a  frees ing  process,  a  semi-infinite 
medium,  initially  in  the  thawed  state  at  70.  suddenly 
undergoes  a  step  change  in  surface  temperature  to  Ts 
which  is  less  than  the  fusion  temperature  T,.  The 
solution  for  the  phase  change  depth  X  is  given  by 

X  =  lyy/t^T  .  (1) 

The  phase  change  parameter  7  can  be  obtained  from 
the  following  equation: 

e1  _  ^  (V  M  be  .  L  yx/iT  (2) 

erf y  Ar,  ( TrTs )  erfeby  C{(TrTs) 

where 

erf>  -  -2- 

V* 

erfey  =  I  -  erf y. 

Carslaw  and  Jaeger  (1959)  present  solutions  to  the 
transcendental  equation  for  7  in  the  case  of  pure  water 
and  ice.  Equation  1 ,  the  Neumann  equation,  is  used 
as  the  basis  for  many  phase  change  studies  and  is  ap¬ 
plicable  to  engineering  problems.  Berggren  (1943)  was 
apparently  the  first  to  actually  apply  the  Neumann  so¬ 
lution  to  soil  phase  change  problems.  Aldrich  and  Payn- 
ter  (1953)  later  used  the  Stefan  form  of  the  phase 
change  solution  to  arrive  at  the  modified  Berggren 
equation.  Equation  1  can  be  changed  to  the  Stefan 
form  as 


X=\yJ{2kv’L)(T{-T%)t  (3) 

where  X,  which  replaces  7,  can  be  determined  from  the 
exact  solution.  Stefan  (1891)  originally  solved  a  simi¬ 
lar  problem  for  the  growth  of  sea  ice  when  the  sensible 
heat  to  latent  heat  ratio  was  small  and  the  water  was 
at  the  freezing  temperature.  Equation  3  reduces  to  the 
Stefan  equation  when  X  equals  one;  hence  the  name 
"Stefan  form." 

The  surface  temperature  of  a  soil  system  does  not 
normally  remain  constant  during  the  freeze  season  and 
the  surface  index  /s  is  often  used. 

X  =  Xx/{2  k,/L]/t.  (4) 

The  surface  index  is  defined  as 
o 

's=  [\TrTf\t,)\dt’  =  [rf-7je  (5) 


where  0  is  the  length  of  the  freeze  season. 

Thus  an  average  constant  surface  temperature  T% 
can  be  calculated  for  the  season  to  be  used  in  cq  1  or 
eq  3,  if  lt  is  known.  Unfortunately  the  surface  index 
is  rarely  available  for  a  location;  however,  the  air  tem¬ 
perature  index  /,  is  usually  tabulated  and  /5  can  be  re¬ 
placed  with  the  n-faclor,  n  defined  as 

n  =  IJIf  or  IJIt  .  (6) 

The  quantity  n  is  the  relation  between  the  air  index 
and  the  surface  index.  A  procedure  for  obtaining  a 
value  at  a  given  site  is  given  by  Lunardini  (1978). 

Finally,  the  modified  Berggren  equation  is  written  as 

X  =  \  .  (7) 

Berg  and  Aitkcn  (1973),  among  many  others,  have 
shown  that  the  modified  Berggren  equation  (eq  7)  gives 
good  results  for  seasonal  phase  change  depths  even  if 
the  surface  temperature  varies  with  time.  The  coefficient 
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y  can  be  found  by  equating  eq  1  and  3  and  substitut¬ 
ing  into  eq  2.  This  leads  to  the  following  equation 
for  y: 

pae^  Xy^  (8) 

erf(X\4T)  erfc(rX\47)  2 

The  parameters  a  and  p  take  into  account  the  soil  tem¬ 
peratures,  specific  heat  and  latent  heat.  The  parameter 
p  is  one-half  of  the  Stefan  number  which  is  the  ratio 
of  the  sensible  heat  and  latent  heat  for  a  soil  system. 
For  small  values  of  p  and  a,  it  can  be  expected  that  X 
will  be  nearly  one  and  eq  7  will  reduce  to  the  Stefan 
equation. 

The  quantities/),  q  and  r  are  ratios  of  the  thermal 
properties  of  the  soil  system  for  the  frozen  and  thawed 


P  =  (V*f)VKf  /«t 
<7  =  Kf/*l 
r  =  V«f/Kt  . 

These  relations  arc  all  for  the  freezing  case.  Aldrich 
and  Payntcr  (1953)  used  the  relations 

a  =  (r0-r,)/(rf-Ts) 

M  =  «rt/z.)(rf-Ts> 

and  noted  that  calculations  with  typical  soil  proper¬ 
ties  indicated  k,/x,  r  1.0,C,/Ct  i  1 ,  and  thus kt/kf 
1 .0.  They  then  solved  eq  8  with p  =  q-r  - 1  and 
obtained  a  widely  used  graph  for  X  (sec  Sanger  1969). 
Actually,  this  procedure  is  only  valid  when  the  water 
content  of  a  soil  is  zero.  Nixon  and  McRobcrts  (1973) 
made  a  parametric  study  of  eq  2,  but  presented  a  graph 
of  X  valid  only  for  r  =  \Jk f,  =  1 .43,  which  was  said  to 
represent  most  soil  conditions. 


SOIL  THERMAL  PROPERTIES 

It  is  clear  that/),  q  and  r  will  vary  for  different  soil 
systems,  but  a  relatively  simple  procedure  can  be  used 
to  generate  these  functions  for  any  soil.  Gold  and 
Lachcnbruch  (1973)  noted  that  the  weighted,  geo¬ 
metric  mean  for  the  thermal  conductivity  of  a  mixture 
gives  results  that  arc  often  as  good  as  more  complicated 
methods.  The  thermal  conductivity  of  a  soil  can  then 
be  expressed  as 

<9> 


where  As,  kv  and  k%  are  the  thermal  conductivities  of 
the  solid,  liquid,  and  gaseous  phases;  *s,  x*,  and  xg  are 
the  volume  fractions  of  the  solid,  liquid,  and  gaseous 
phases.  For  soil  systems,  the  thermal  conductivity  of 
the  solids  and  gases  will  not  vary  significantly  as  phase 
change  occurs  (Kerstcn  1949).  There  will  be  only  a 
small  error  if  it  is  assumed  that  the  frozen  state  con¬ 
tains  only  ice  with  no  unfrozen  water.  Thus  the  ratio 
of  the  frozen  to  unfrozen  conductivities  of  the  soil 
mixture  can  be  related  to  the  thermal  conductivity  of 
ice  and  water  as  follows: 

V*i  =  00) 

where &w,  hi  are  the  thermal  conductivities  of  water 
and  ice.  The  volumetric  specific  heat  for  the  system 
may  be  expressed  for  the  thawed  and  frozen  states  as 
follows: 

Ct  =Csi(1-*v.)+Cw*t  (11) 

Cf  =Cft(1-xf)  +  Cix,  (12) 

where  Csl,  Cs)  =  specific  heats  of  unfrozen  and  frozen 
soil  solids.  The  neglect  of  the  gas  phase  is  insignificant 
since  the  density  of  the  gas  (air)  is  low.  It  is  fortunate 
that  the  specific  heats  of  different  soil  solids  and  ice 
are  all  similar  in  magnitude.  For  example,  the  volu¬ 
metric  specific  heal  of  organic  solids  is  about  2300  k ) / 
m*  K,  for  mineral  solids  it  is  1760,  while  for  ice  it  is 
1920  (sec  Lunardini  1971).  The  properties  of  the  fro 
zen  materials  are  evaluated  at  25°F  (269  K)  while  the 
thawed  values  arc  at  40°K  (277.4  K).  If  one  assumes 
that  the  values  for  the  solids,  except  for  ice,  change 
little  through  the  phase  change  then 

C,/C,  =  1  +  |(Cw/C,)-1|av  (13) 


C,/Ct  =  1  +  1.023  v 

These  relations  indicate  that  the  soil  property  ratios 
may  vary  as 

0.25  ■  A-,/*,  .  1.0 

1.0  C,/C,  <  2.1 

0.13  k,/k,  .  1.0. 

The  property  values  to  use  in  eq  8  can  then  be  expressed 
as  simple  functions  of  the  soil  water  content 

q  -  Rv  (1*  1.023**) 
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Figure  3.  Freezing  case,  Neumann  equation  parameter,  =  0.8. 


Figure  4.  Freezing  case,  Neumann  equation  parameter,  xv  -  1.0. 


r=  y/q 
P  =  r/RXt. 

The  ratio  of  the  thermal  conductivity  of  ice  to 
water  R  is  about  3.89.  Equation  8  can  be  solved  nu¬ 
merically  to  find  the  roots,  which  are  the  values  of  X. 
Figures  1  through  4  give  the  values  of  X  to  use  in  eq  3 
or  eq  7  for  the  freezing  case. 

THAW  CASE 

At  first  glance,  it  might  appear  that  the  same  rela¬ 
tions  could  be  used  for  either  the  thawing  or  the  freez¬ 
ing  case.  This,  however,  is  not  true.  In  the  thawing 
case,  the  medium  is  initially  frozen  at  Ta  and  energy 
must  be  conducted  through  the  thawed  layer  from  the 
phase  change  interface.  Since  the  thermal  conductivity 
of  the  thawed  region  is  considerably  less  than  that  of 
the  frozen,  the  heat  flow  will  be  reduced,  even  with 
the  same  temperature  gradient.  However,  the  general 
form  of  the  equation  will  be  the  same;  after  making 
appropriate  changes  for  the  property  values.  The  thaw 
depth  is  expressed  as 

X  -  \yj(2kt  lxn)/L  (14) 


where  X  is  again  given  by  eq  8  but  the  parameters  a 
and  p  are  now  a  =  (Tf-T0)l(  7 \-Tt)  and  p  =  (CJ2L) 

( rs -  r, ) .  The  functionsp,  q  andr  all  change  because 
of  the  property  changes  of  the  thawed  and  frozen  states. 

q  =  1/[/?'l(l  +  1 .023  ) | . 

r=\Jq  ■ 

p  =  rRe. 

The  X  values  for  thawing  arc  now  given  in  Figures  5 
through  7.  Notice  that  when  =  0  the  X  values  are 
the  same  for  freezing  or  thawing  (Fig.  1 ),  which  is  the 
Aldrich  (1953)  case  and  the  Sanger  (1969)  graph.  The 
charts  for  X  are  for  the  exact  solution  of  the  Neumann 
problem  with  property  values  typical  of  soil  systems. 

DISCUSSION 

A  comparison  of  the  value  of  X  obtained  with  Fig¬ 
ures  1  through  1  and  the  value  obtained  from  the  chart 
given  in  Sanger  (1969)  shows  that  the  values  can  differ 
by  at  least  ±10%.  Also,  for  the  same  values  of  xt,  a 
and  p,  the  value  of  X  differs  by  ±  10%  when  the  freez¬ 
ing  and  thawing  cases  arc  considered.  Since  Sanger's 
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Figure  5.  Thawing  ease,  Neumann  equation  parameter  xv  ■  0.4. 
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chart  uses  property  ratios  of  one,  there  is  no  distinction 
between  thawing  and  freezing.  While  these  variations 
are  not  extreme  they  will  lead  to  the  same  percentage 
differences  for  the  calculated  depths  of  freeze  or  thaw. 
No  additional  computational  work  is  involved  and  it 
seems  desirable  to  use  the  graphs  presented  here. 

If  the  modified  Berggrcn  equation  is  programmed 
into  a  computer  then  eq  8  would  be  used  directly  to 
generate  the  X  values.  Equation  8  can  be  solved  by 
any  convenient  scheme  such  as  iteration  or  Newton’s 
method.  The  values  for  Figures  1  through  7  were  ob¬ 
tained  with  Newton’s  method.  The  quantities p,  q  and 
r  can  be  calculated  directly  if  the  thermal  property 
ratios  are  specified  or  by  using  the  functional  depend¬ 
ence  upon  the  volumetric  water  content  derived  in  this 
report. 
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